% X = [X1 X2 ... Xn]

function Y = fastCov(X)
X_bar = mean(X,2);
X1=bsxfun(@minus,X,X_bar);
Y=X1*X1'/(size(X,2)-1);